﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using RayDen.Library.Core;
using RayDen.Library.Core.Primitives;

namespace Lt.Engine.Types
{
    public class RayIntersection
    {
        public Point Hitpoint;
        public Normal GeoNormal, ShadingNormal;
        public Vector dpdu, dpdv, dpdx, dpdy;
        public UV TextureCoords, BaryCentricPoint;

        public float dudx, dvdx, dudy, dvdy;

        public void ComputeDifferential(RayDifferential ray)
        {

            if (ray.hasDifferentials)
            {
                // Estimate screen space change in $\pt{}$ and $(u,v)$

                // Compute auxiliary intersection points with plane
                var p = new Normal(Hitpoint.x, Hitpoint.y, Hitpoint.z);
                float d = -Normal.Dot(ref GeoNormal, ref p);
                var rxv = new Vector(ray.rxOrigin.x, ray.rxOrigin.y, ray.rxOrigin.z);
                float tx = -(Vector.Dot(ref GeoNormal, ref rxv) + d) / Vector.Dot(ref GeoNormal, ref ray.rxDirection);
                if (float.IsNaN(tx))
                {
                    dudx = dvdx = 0f;
                    dudy = dvdy = 0f;
                    dpdx = dpdy = new Vector(0f);
                    return;
                }
                Point px = ray.rxOrigin + tx * ray.rxDirection;
                Vector ryv = new Vector(ref ray.ryOrigin);
                float ty = -(Vector.Dot(ref GeoNormal, ref ryv) + d) / Vector.Dot(ref GeoNormal, ref ray.ryDirection);
                if (float.IsNaN(ty))
                {
                    dudx = dvdx = 0f;
                    dudy = dvdy = 0f;
                    dpdx = dpdy = new Vector(0f);
                    return;
                }

                Point py = ray.ryOrigin + ty * ray.ryDirection;
                dpdx = px - Hitpoint;
                dpdy = py - Hitpoint;

                // Compute $(u,v)$ offsets at auxiliary points

                // Initialize _A_, _Bx_, and _By_ matrices for offset computation
                float[][] A = new float[2][];
                float[] Bx = new float[2];
                float[] By = new float[2];
                int[] axes = new int[2];
                if (Math.Abs(GeoNormal.x) > Math.Abs(GeoNormal.y) && Math.Abs(GeoNormal.x) > Math.Abs(GeoNormal.z))
                {
                    axes[0] = 1; axes[1] = 2;
                }
                else if (Math.Abs(GeoNormal.y) > Math.Abs(GeoNormal.z))
                {
                    axes[0] = 0; axes[1] = 2;
                }
                else
                {
                    axes[0] = 0; axes[1] = 1;
                }

                // Initialize matrices for chosen projection plane
                A[0][0] = dpdu[axes[0]];
                A[0][1] = dpdv[axes[0]];
                A[1][0] = dpdu[axes[1]];
                A[1][1] = dpdv[axes[1]];
                Bx[0] = px[axes[0]] - p[axes[0]];
                Bx[1] = px[axes[1]] - p[axes[1]];
                By[0] = py[axes[0]] - p[axes[0]];
                By[1] = py[axes[1]] - p[axes[1]];
                if (!MathLab.SolveLinearSystem2X2(A, Bx, out dudx, out dvdx))
                {
                    dudx = 0f; dvdx = 0f;
                }
                if (!MathLab.SolveLinearSystem2X2(A, By, out dudy, out dvdy))
                {
                    dudy = 0f; dvdy = 0f;
                }
            }
            else
            {
                dudx = dvdx = 0f;
                dudy = dvdy = 0f;
                dpdx = dpdy = new Vector(0f);
            }
        }
    }
}
